#include "rsdecode.h"

RSdecode::RSdecode()
{


    RS_PERM_ERROR = -1;
    RS_CORRECT_ERROR = -2;

    POLYNOMIAL = 0x1d;
}

void RSdecode::setNpar(int a)
{
    npar = a;
}

int RSdecode::calcSigmaMBM(int *sigma, int *omega, int *syn)
{
          int sg0[npar];
          int sg1[npar];
          sg0[1] = 1;
          sg1[0] = 1;
            int jisu0 = 1;
            int jisu1 = 0;
            int m = -1;

            for(int n = 0; n < npar; n++) {
                // 判別式を計算
                int d = syn[n];
                for(int i = 1; i <= jisu1; i++) {
                    d ^= mul(sg1[i], syn[n - i]);
                }
                if(d != 0) {
                    int logd = toLog(d);
                    int wk[npar];
                    for(int i = 0; i <= n; i++) {
                        wk[i] = sg1[i] ^ mulExp(sg0[i], logd);
                    }
                    int js = n - m;
                    if(js > jisu1) {
                        m = n - jisu1;
                        jisu1 = js;
                        if(jisu1 > npar / 2) {
                            return -1;				// σの次数がnpar / 2を超えたらエラー
                        }
                        for(int i = 0; i <= jisu0; i++) {
                            sg0[i] = divExp(sg1[i], logd);
                        }
                        jisu0 = jisu1;
                    }

                    for (int c = 0; c < sizeof (wk); c++)
                        sg1[c] = wk[c];
                }
               // System.arraycopy(sg0, 0, sg0, 1, Math.min(sg0.length - 1, jisu0));
                int min = (sizeof(sg0) - 1 < jisu0)? sizeof(sg0) - 1 : jisu0;
                int j = 0;
                for (int i = 1; i < min; i++) {
                    sg0[j++] = sg0[i];
                }

                sg0[0] = 0;
                jisu0++;
            }
            mulPoly(omega, sg1, syn);
            //System.arraycopy(sg1, 0, sigma, 0, Math.min(sg1.length, sigma.length));
            int min2 = (sizeof(sg1) - 1 < sizeof(sigma))? sizeof(sg1) - 1 : sizeof(sigma);
            int k = 0;
            for (int i = 1; i < min2; i++) {
                sigma[k++] = sg1[i];
            }
            return jisu1;
}

int RSdecode::chienSearch(int *pos, int n, int jisu, int *sigma)
{
    /* * σ(z) = (1-α^i*z)(1-α^j*z)(1-α^k*z)
         *       = 1 + σ1z + σ2z^2 +...
         * σ1 = α^i + α^j + α^k
         * 上記の性質を利用して、プチ最適化
         * last = σ1から、見つけた解を次々と引いていくことにより、最後の解はlastとなる
         */
        int last = sigma[1];

        if(jisu == 1) {
            // 次数が1ならば、lastがその解である
            if(toLog(last) >= n) {
                return RS_CORRECT_ERROR;	// 範囲外なのでエラー
            }
            pos[0] = last;
            return 0;
        }

        int posIdx = jisu - 1;		// 誤り位置格納用インデックス
        for(int i = 0; i < n; i++) {
            /*
             * σ(z)の計算
             * w を1(0乗の項)に初期化した後、残りの項<1..jisu>を加算
             * z = 1/α^i = α^Iとすると
             * σ(z) = 1 + σ1α^I + σ2(α^I)^2 + σ3(α^I)^3 + ... + σ<jisu>/(α^I)^<jisu>
             *       = 1 + σ1α^I + σ2α^(I*2) + σ3α^(I*3) + ... + σ<jisu>α^(I*<jisu>)
             */
            int z = 255 - i;					// z = 1/α^iのスカラー
            int wk = 1;
            for(int j = 1; j <= jisu; j++) {
                wk ^= mulExp(sigma[j], (z * j) % 255);
            }
            if(wk == 0) {
                int pv = toExp(i);		// σ(z) = 0の解
                last ^=  pv;					// lastから今見つかった解を引く
                pos[posIdx--] = pv;
                if(posIdx == 0) {
                    // 残りが一つならば、lastがその解である
                    if(toLog(last) >= n) {
                        return RS_CORRECT_ERROR;	// 最後の解が範囲外なのでエラー
                    }
                    pos[0] = last;
                    return 0;
                }
            }
        }
        // 探索によりデータ長以内に、jisu個の解が見つからなかった
        return RS_CORRECT_ERROR;
}

void RSdecode::doForney(int *data, int length, int jisu, int *pos, int *sigma, int *omega)
{
    for(int i = 0; i < jisu; i++) {
                int ps = pos[i];
                int zlog = 255 - toLog(ps);					// zのスカラー

                // ω(z)の計算
                int ov = omega[0];
                for(int j = 1; j < jisu; j++) {
                    ov ^= mulExp(omega[j], (zlog * j) % 255);		// ov += ωi * z^j
                }

                // σ'(z)の値を計算(σ(z)の形式的微分)
                int dv = sigma[1];
                for(int j = 2; j < jisu; j += 2) {
                    dv ^= mulExp(sigma[j + 1], (zlog * j) % 255);	// dv += σ<j+1> * z^j
                }

                /*
                 * 誤り訂正 E^i = α^i * ω(z) / σ'(z)
                 * 誤り位置の範囲はチェン探索のときに保証されているので、
                 * ここではチェックしない
                 */
                data[toPos(length, ps)] ^= mul(ps, div(ov, dv));
            }
}

int RSdecode::decode(int *data, int length, bool noCorrect)
{
    if(length < npar || length > 255) {
                return RS_PERM_ERROR;
            }
            // シンドロームを計算
            int syn[npar];
            if(calcSyndrome(data, length, syn)) {
                return 0;		// エラー無し
            }
            // シンドロームよりσとωを求める
            int sigma[npar / 2 + 2];
            int omega[npar / 2 + 1];
            int jisu = calcSigmaMBM(sigma, omega, syn);
            if(jisu <= 0) {
                return RS_CORRECT_ERROR;
            }
            // チェン探索により誤り位置を求める
            int pos[jisu];
            int r = chienSearch(pos, length, jisu, sigma);
            if(r < 0) {
                return r;
            }
            if(!noCorrect) {
                // 誤り訂正
                doForney(data, length, jisu, pos, sigma, omega);
            }
            return jisu;
}

int RSdecode::decode(int *data, int length)
{
    return decode(data, length, false);
}

int RSdecode::decode(int *data)
{
    return decode(data, sizeof (data), false);
}

void RSdecode::initGaloisTable()
{
    int d = 1;
    for(int i = 0; i < 255; i++) {
                expTbl[i] = expTbl[255 + i] = d;
                logTbl[d] = i;
                d <<= 1;
                if((d & 0x100) != 0) {
                    d = (d ^ POLYNOMIAL) & 0xff;
                }
    }
}

void RSdecode::mulPoly(int *seki, int *a, int *b)
{
    for (int i = 0 ; i < sizeof (seki); i++) {
        seki[i] = 0;
    }

    for(int ia = 0; ia < sizeof (a); ia++) {
                if(a[ia] != 0) {
                    int loga = logTbl[a[ia]];
                    int ib2 = ( sizeof (b) < sizeof (seki) - ia)? sizeof (b) : sizeof (seki) - ia ;
                    for(int ib = 0; ib < ib2; ib++) {
                        if(b[ib] != 0) {
                            seki[ia + ib] ^= expTbl[loga + logTbl[b[ib]]];	// = a[ia] * b[ib]
                        }
                    }
                }
     }

}

bool RSdecode::calcSyndrome(int *data, int length, int *syn)
{
    int hasErr = 0;
    for(int i = 0; i < sizeof (syn);  i++) {
              int wk = 0;
              for(int idx = 0; idx < length; idx++) {
                     wk = data[idx] ^ ((wk == 0)? 0 : expTbl[logTbl[wk] + i]);		// wk = data + wk * α^i
              }
              syn[i] = wk;
              hasErr |= wk;
    }

    return hasErr == 0;
}

int RSdecode::toExp(int a)
{
    return expTbl[a];
}

int RSdecode::toLog(int a)
{
    return logTbl[a];
}

int RSdecode::toPos(int lenght, int a)
{
    return lenght - 1 - logTbl[a];
}

int RSdecode::mul(int a, int b)
{
    return (a == 0 || b == 0)? 0 : expTbl[logTbl[a] + logTbl[b]];
}

int RSdecode::mulExp(int a, int b)
{
    return (a == 0)? 0 : expTbl[logTbl[a] + b];
}

int RSdecode::div(int a, int b)
{
    return (a == 0)? 0 : expTbl[logTbl[a] - logTbl[b] + 255];
}

int RSdecode::divExp(int a, int b)
{
    return (a == 0)? 0 : expTbl[logTbl[a] - b + 255];
}

int RSdecode::inv(int a)
{
    return expTbl[255 - logTbl[a]];
}

